
COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
   lpstr,rho_hi,rho_h2,rho_star


restore,'./irdc_dist_model/bgps_rb3.sav',/ver

;; Do some work, and place this info into the structure!
n       = n_elements(rb3)
tau     = dblarr(n)
f_data  = dblarr(n)
sig_f   = dblarr(n)
nearfar = bytarr(n)

myps,'./irdc_dist_model/analysis_plots/compare_sigma_cont.eps',$
     xsize=10,ysize=13.3

sym_tau = cgSymbol('tau')
sym_sig = cgSymbol('sigma')
sym_gam = cgSymbol('Gamma')

multiplot,[4,4],ygap=0.04,/doxaxis


;; Loop though cases
FOR cas = 1, 4 DO BEGIN
   
   cstr = string(cas,format="(I0)")
   
   restore,'./irdc_dist_model/data/case'+cstr+$
           '_irdc_TdFIT.sav'
   
   Td = res[4]
   TEMP_DEPENDENT_VARIABLES, Td
   gamma = FLUX2TAU_BGPS(Td, /mJy, /apper)
   
   imin = rb3.i_min
   imir = rb3.i_mir
   smin = rb3.sig_imin
   smir = rb3.sig_imir
   sess = rb3.s_err
   C = rb3.c_meas
   
   var_min = (imin*imin/imir/imir)*(smin*smin/imin/imin)
   var_mir = (imin*imin/imir/imir)*(smir*smir/imir/imir)
   var_ess = sess*sess*exp(-2.d*tau)*C*C*gamma*gamma/$
             ((1.d - exp(-tau))*(1.d - exp(-tau)))
   
   sum = var_min + var_mir + var_ess
   
   f_min = var_min / sum
   f_mir = var_mir / sum
   f_ess = var_ess / sum
   f_sum = sum / sum
   
   plotsym,0,0.25,/fill
   xlog = 1b
   IF cas EQ 1 THEN ytit='Fraction of total' ELSE ytit=''
   IF cas EQ 3 THEN xtit='Uncertainty in f!ddata!n' ELSE xtit=''
   
   plot,sig_f,f_ess,/nodata,xr=[0.011,10.],/xst,yr=[-0.02,1.02],/yst,xlog=xlog,$
        ytit=ytit,xtit=xtit
   oplot,sig_f,f_ess,psym=8,color=cgColor('Dodger Blue')
   oplot,sig_f,f_min,psym=8,color=cgColor('Red')
   oplot,sig_f,f_mir,psym=8,color=cgColor('Dark Green')
   
   ;;linfit_plot,sig_f,f_ess,xlog=xlog,/legend,/top,/left
   
   
   al_legend,/top,/right,box=0,textcolor='Black',charsize=0.8,$
             ['Case '+cstr,'T!dd!n = '+string(Td,format="(F0.1)")+'K']
   
   multiplot,/doxaxis
   
ENDFOR

;; Middle row of figure...

imin = rb3.i_min
imir = rb3.i_mir
smin = rb3.sig_imin
smir = rb3.sig_imir
sess = rb3.s_err
C = rb3.c_meas

plotsym,0,0.25,/fill

ytit=sym_sig+'!dImin!n / I!dMIR!n' 
xtit='I!dmin!n [MJy sr!u-1!n]'

y = smin/imir
sc_min = y

plot,imin,y,psym=8,ytit=ytit,xtit=xtit,/ylog,/xlog,$
     ytickformat='exponent10',/nodata
oplot,imin,y,psym=8,color=cgColor('Red')
linfit_plot,imin,y,/legend,/top,/left,charsize=0.4,/ylog,/xlog,$
            color=cgColor('Red')

vline,0.1,/horiz,color=cgColor('red6'),thick=5,linestyle=3,/log

multiplot

multiplot,/doxaxis,/doyaxis

ytit='I!dmin!n '+sym_sig+'!dIMIR!n / I!dMIR!u2!n' 
xtit='I!dMIR!n [MJy sr!u-1!n]'

y = imin*smir/imir/imir
sc_mir = y

plot,imir,y,psym=8,ytit=ytit,xtit=xtit,/ylog,/xlog,$
     ytickformat='exponent10',/nodata
oplot,imir,y,psym=8,color=cgColor('Dark Green')
linfit_plot,imir,y,/legend,/top,/left,charsize=0.4,/ylog,/xlog,$
            color=cgColor('Dark Green')

vline,0.1,/horiz,color=cgColor('red6'),thick=5,linestyle=3,/log

multiplot

multiplot,/doxaxis,/doyaxis

;; Next row...


;; Loop though cases
FOR cas = 1, 4 DO BEGIN
   
   cstr = string(cas,format="(I0)")
   
   restore,'./irdc_dist_model/data/case'+cstr+$
           '_irdc_TdFIT.sav'
   
   Td = res[4]
   TEMP_DEPENDENT_VARIABLES, Td
   gamma = FLUX2TAU_BGPS(Td, /mJy, /apper)
   
   imin = rb3.i_min
   imir = rb3.i_mir
   smin = rb3.sig_imin
   smir = rb3.sig_imir
   sess = rb3.s_err
   s_40 = rb3.s_peak
   C = rb3.c_meas
   tau = gamma * rb3.s_peak
   
   sig_tau = gamma*sess
   

   sum = var_min + var_mir + var_ess
   
   plotsym,0,0.25,/fill
   xlog = 1b
   IF cas EQ 1 THEN $
      ytit='e!u-'+sym_tau+'!nC'+sym_gam+sym_sig+'!dS!n / (1 - e!u-'+$
           sym_tau+'!n)' ELSE ytit=''
   IF cas EQ 3 THEN xtit='S [mJy]' ELSE xtit=''
   
   y = exp(-tau)*C*gamma*sess/(1.d - exp(-tau))
   sc_ess = y
   ;;print,m4_stat(s_40)
   
   plot,s_40,y,/nodata,/ylog,yr=[1.d-5,0.5],xr=[50,1.d5],$
        /xst,/yst,xlog=xlog,ytit=ytit,xtit=xtit
   oplot,s_40,y,psym=8,color=cgColor('Dodger Blue')
   
   linfit_plot,s_40,y,/xlog,/ylog,/legend,/top,/left,$
               color=cgColor('Dodger Blue'),charsize=0.5
   
   
   al_legend,/top,/right,box=0,textcolor='Black',charsize=0.8,$
             ['Case '+cstr,'T!dd!n = '+string(Td,format="(F0.1)")+'K']
   
   vline,0.1,/horiz,color=cgColor('red6'),thick=5,linestyle=3,/log
   
   multiplot,/doxaxis
   
ENDFOR

;;Bottom row

plot,sc_mir,sc_min,psym=8,/xlog,/ylog,ytickformat='exponent10',$
     xtickformat='exponent10',ytit=sym_sig+'!dImin!n / I!dMIR!n',$
     xtit='I!dmin!n '+sym_sig+'!dIMIR!n / I!dMIR!u2!n',xr=[1d-5,0.999],/xst
linfit_plot,sc_mir,sc_min,/xlog,/ylog,/legend,/top,/left,charsize=0.5
plots,[1d-3,1d-3]
plots,[1,1],/cont,linestyle=3,thick=3,color=cgColor('TG3')

multiplot,/doxaxis

plot,sc_ess,sc_min,psym=8,/xlog,/ylog,xtickformat='exponent10',$
     xtit='e!u-'+sym_tau+'!nC'+sym_gam+sym_sig+'!dS!n / (1 - e!u-'+$
           sym_tau+'!n)'
linfit_plot,sc_ess,sc_min,/xlog,/ylog,/legend,/top,/left,charsize=0.5
plots,[1d-3,1d-3]
plots,[1,1],/cont,linestyle=3,thick=3,color=cgColor('TG3')

multiplot

multiplot,/doxaxis,/doyaxis

plot,sc_ess,sc_mir,psym=8,/xlog,/ylog,ytickformat='exponent10',$
     xtickformat='exponent10',ytit='I!dmin!n '+sym_sig+$
     '!dIMIR!n / I!dMIR!u2!n',xtit='e!u-'+sym_tau+'!nC'+sym_gam+sym_sig+$
     '!dS!n / (1 - e!u-'+sym_tau+'!n)'
linfit_plot,sc_ess,sc_mir,/xlog,/ylog,/legend,/top,/left,charsize=0.5
plots,[1d-5,1d-5]
plots,[1,1],/cont,linestyle=3,thick=3,color=cgColor('TG3')






multiplot,/reset
multiplot,/default
multiplot,/reset

myps,/done

END
